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Abstract. We examine the out-of-cquilibrium dynamical evolution of density profiles 
of ultrasoft particles under time-varying external confining potentials in three spatial 
dimensions. The theoretical formalism employed is the dynamical density functional 
theory (DDFT) of Marini Bettolo Marconi and Tarazona [J. Chem. Phys. 110, 8032 
(1999)], supplied by an equilibrium excess free energy functional that is essentially 
exact. We complement our theoretical analysis by carrying out extensive Brownian 
Dynamics simulations. We find excellent agreement between theory and simulations 
for the whole time evolution of density profiles, demonstrating thereby the validity of 
the DDFT when an accurate equilibrium free energy functional is employed. 

Density functional theory (DFT) is a very powerful tool for the quantitative 
description of the equilibrium states of many-body systems under arbitrary external 
fields. It rests on the exact statement that the Helmholtz free energy of the system, 
F[p], is a unique functional of the inhomogeneous one-particle density p(r). Moreover, 
the equilibrium profile Po( r ) minimises F[p] under the constraint of fixed particle number 
N p. The task is then to approximate the unknown functional F[p] from which all 
equilibrium properties of the system follow. Much more challenging is the problem 
of studying out- of equilibrium dynamics of many-body systems, for which analogous 
uniqueness and minimisation principles are lacking. In this paper, we present results 
based on a recently-proposed dynamical density functional theory (DDFT) formalism 
and we demonstrate that the latter is capable of describing out-of-equilibrium diffusive 
processes in colloidal systems at the Brownian time scale. 

We are concerned with the dynamics of typical soft-matter systems, such as 
suspensions of mesoscopic spheres and polymer chains in a microscopic solvent {2J . The 
enormous difference in the masses of the suspended particles and the solvent molecules 
implies a corresponding separation in the relaxational time scales of the two. At times 
of the order of the Fokker-Planck scale, rpp ~ 10 -14 sec, the solvent coordinates are long 
relaxed to thermal equilibrium. On the Brownian diffusive time scale, tb ~ 10 -9 sec, 
the momentum coordinates of the solute particles relax to equilibrium with the heat 
bath of the solvent molecules and thus a statistical description involving only the 
positions of the colloids is feasible >3]. In this regime, the evolution of the coordinates 
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{r%(t), r 2 (t), . . . , r N (t)} of the N colloidal particles is described by the set of stochastic 
Langevin equations: 



dr*(i) 



dt 



-rv r 



Wi(t). 



(1) 



£v(|r,-r,,|)+Ve*(r<,t) 

In Eq. above, V(\vi — Yj\) is the pair (effective) interaction potential between 
the mesoscopic particles 0, V r ext (rj,t) is the externally acting potential and Wj(t) = 
[wf (t), w\ (£), w-(t)} is a stochastic term representing the random collisions with the 
solvent molecules and having the properties: 



«(*))= and {w < *{t)w p j {t')) = 2D8 ij 5 a p8{t-t l ), (2) 

where the averages (• • •) are over the Gaussian noise distribution and a, (3 — x, y, z, 
the Cartesian components. The constants T and D stand for the mobility and 
diffusion coefficients of the particles, respectively, and the Einstein relation gives 
T/D = (/cgT) -1 = (3. Applying the rules of the Ito stochastic calculus, Marconi and 
Tarazona jU IHj recasted the above equations into the form 

x dp{r,t) _ , 



dt 



+ V r 



k B TV rP (r, t) + p(r, t)V r V ext (r, t) 
dV (p(r,t)p(r',t))V r F(r-r / ; 



(3) 



Here, p(r,t) = J2i$(. r i(t) ~~ r ) i s the usual one-particle density operator and p(r,t) = 
(p(r,t)) is the noise-average of this quantity. Up to this point, all is exact. Now, the 
following physical assumption (A) is introduced: as the system follows its relaxative 
dynamics, the instantaneous two-particle correlations can be approximated with those 
of a system in thermodynamic equilibrium with the same, static one-particle density 
p(r) as the noise-averaged dynamical one-particle density p(r,t). Then, Eq. (jHJ) can be 
cast into a form involving exclusively the equilibrium density functional F[p] as [HE] 



(4) 



with the free energy functional 

F[p}=k B T J d 3 rp(r) {in [p(r)A 3 ] -l}+F cx [p] + J d 3 rV ext (r, t) p(r). (5) 

The dynamical equation of motion (j3J) was in fact first derived in a phenomenological 
way by Dieterich, Frisch and Majhofer jB]. 

In carrying out concrete calculations with the theory put forward above and 
in comparing them with Brownian Dynamics (BD) simulation results based on the 
microscopic equations of motion, Eq. (JTJ), two sources of possible discrepancies exist: 
first, the fundamental assumption (A) and second the approximate nature of the 
equilibrium density functional F ex [p] of Eq. (jSJ). In this work we focus our attention to 
ultrasoft particles for which a very accurate and simple functional F[p] is known, namely 
the mean-field or random-phase approximation (RPA) functional given by Eq. (JHJ) below. 
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This guarantees that one can explore the accuracy of the fundamental assumption (A) 
under well-defined external conditions. 

Consider a one-component system of ultrasoft particles. It has been demonstrated 
that for such systems the following RPA-functional is quasi-exact [71 IHj 1§1 ITTH ITT) ITS } ITS] : 

F cx [p] = \J J d 3 rdW(|r - r'\)p(v)p{v'). (6) 

Eq. (J3]) takes now with the help of Eqs. (0) and © the form 

r -i5eM = fcfl r^(r,t) + V r p(r,t). J dVV r V(|r-r'|)p(r',t) 

+ p(r,t) J dVV^(|r-r'|)p(r',t) 

+ V r p(r, t) ■ V r Kxt(r, t) + p(r, t)V 2 r V ext (r, t). (7) 

Given an initial density field p(r, t = 0) and a prescribed external potential Kxt( r ) £), Eq. 
(J7J can be numerically solved to yield p(r, t). In this work we apply an ultrasoft Gaussian 
pair potential between the interacting particles that has been shown to describe the 
effective interaction between the centres of mass of polymer chains in athermal solvents 

V(r) = eexp[-(r/a) 2 ]. (8) 

We set e = ksT providing the energy unit for the problem, whereas a, which corresponds 
to the gyration radius of the polymers, will be the unit of length henceforth. Accordingly, 
the natural time scale of the problem, providing the unit of time in this work, is 
the Brownian time scale tb = a 2 /{eT). Eq. (jJJ is solved using standard numerical 
techniques, and for a variety of time-dependent confining external potentials V ex t(r,t) 
to be specified below. 

Brownian Dynamics simulations of Eq. (Q) are also straightforward to carry out. 
The Langevin equations of motion including the external field are numerically solved 
using a finite time-step At = 0.003 tb in all simulations, and the technique of Ermak 
^HlQSj- In order to obtain the time-dependent density p(r, t) we perform a large number 
N run of independent runs, typically N Tun = 5000, and average the density profile over 
all configurations for a fixed time t. 

We focus to external fields that correspond to a sudden change, i.e., V ext (r,t) = 
$i(r)0(— t) + $ 2 ( r )@(^)- These force the system to relax from the equilibrium density 
pi(r) = p(r, t < 0), compatible to the external potential $i(r), to the new equilibrium 
density P2(r) = p(r, t — > oo), corresponding to the external potential $2( r )- Important 
questions related to such processes are what is the typical relaxation time r for such a 
procedure and how does the system cross over from one equilibrium density to the other. 
We consider two kinds of confinements: spherical ones, V ext (r,t) = V ext (r, t), where 
r = |r|, and planar ones between two walls perpendicular to the ^-Cartesian direction, 
Kxt( r , t) = V ext (z,t). In these cases we obtain p(r, t) = p(r,t) and p(r, t) = p(z,t), 
correspondingly, and the solution of Eq. (J7J is greatly simplified since the integrals take 
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Figure 1. DDFT (solid lines) and BD (noisy lines) results for the time development 
of the radial density profiles p(r) in the spherical confining potential V^t ( r ' with (a) 
Ri = 4.0 and i?2 = 6.0 and (b) i?i = 6.0 and R2 = 4.0. The shown profiles are for 
the times to = 0, t\ = 0.06, t 2 = 0.18, £3 = 0.54 and £4 = 2.0, all in re-units. The last 
time is practically equivalent to t — 00, since the system there has fully relaxed into 
equilibrium. In all figures, red curves denote the initial and blue ones the final static 
profile. 



the form of one-dimensional convolutions that can be evaluated very rapidly by use of 
fast Fourier transform techniques. 

Three different external confinements have been specifically investigated. Two 
spherical ones 

V%(r,t) = % [(r/R^&i-t) + (r/R 2 y&(t)} ; (9) 
V<S(r,t) = $0 [(r/R 1 ) w e(-t) + (r/R 2 ) w e(t)] ; (10) 
and one slab confinement 

v®(z,t) = % [(z/z 1 ) 10 e(-t) + <p(z/z 2 ) 10 e(t)] . (11) 

The energy scale $0 se ts the strength of the confining potential and is fixed to 
$0 = 10 ksT for all three confinements. The only difference between the external 
potential for times t < and for t > lies in the different length scales R 2 7^ Ri for Eqs. 
(JHJ) and (jlOj) . and Z\ 7^ Z 2 for Eq. (jlljl . For each confinement we consider two cases 
that give rise to two different dynamical processes: R\ < R 2 (Z 1 < Z 2 ), enforcing an 
expansion of the system and R\ > R 2 {Z\ > Z 2 ), bringing about a compression of the 
same. For the spherical confinements an additional parameter is the particle number 
N = j d 3 rp(r,t) which is a conserved quantity, as is clear from Eq. (@J) that has the 
form of a continuity equation. N enters the formalism through the normalisation of the 
density field at t = 0. For both spherical confinements the particle number is N = 100. 
In the slab confinement, Eq. (jlljl . the conserved quantity is the density per unit area 
Po = J dz p(z,t). We choose p^a 2 = 10. In all cases examined, the typical relaxation 
time was found to be of order tb; after typically t = 2tb, the system fully relaxes into 
the new equilibrium profile. 

In Fig. we show the results for the harmonic confining potential of Eq. Q- It can 
be seen that the theory reproduces the time evolution of the density profile, both for the 
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Figure 2. DDFT (solid lines) and BD (noisy lines) results for the time development 
of the radial density profiles p(r) in the spherical confining potential Vg Xt (r, t) with (a) 
Ri = 4.0 and i?2 = 5.0 and (b) R\ — 5.0 and R2 = 4.0. The shown profiles are for the 
times t a = 0, t x = 0.03, t 2 = 0.06, t 3 = 0.12, t A = 0.24, t 5 = 0.48 and t 6 = 2.0 (in units 
of r B ). 



expansion [Fig. ^a)] and the compression [Fig. EJb)] processes. An asymmetry in the 
two processes can be already discerned: the compression is not the 'time reversed' of the 
expansion and this effect will be much stronger in the examples to follow. Though the 
profiles of the system are considerably different from those of an ideal gas, i.e., effects 
of the interparticle interaction are present, the confining potential is smooth enough, so 
that the profiles are devoid of pronounced correlation peaks. 

The situation is different for the external potential of Eq. (jlOj) . Here, the power- 
law dependence is much more steep, so that the Gaussian fluid develops correlation 
peaks close to the 'walls' of the confining field. The dynamical development of the 
profiles for the expansion and compression processes are shown in Fig. El Here, the 
asymmetry between the expansion and the compression processes is evident. In the 
former case, seen in Fig. Efa), the expansion of the confining potential leaves behind 
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Figure 3. The second moment of the radial density profile, m,2(t), defined in Eq. Ijf 2(1 
plotted against the time t for the spherical confinement V^, xt (r, t) . Circles correspond 
to radii R\ = 4.0 and i?2 = 5.0 (expansion) and squares show the resulting curve for 
the inverse process, R\ — 5.0 and R2 — 4.0 (compression). The lines are the analytical 
fits shown in the text. Solid line: Eq. (|f 311 : long-dashed line: Eq. I|14fl . The arrows 
mark the characteristic time scales defined in these two equations. 



a density profile that has very strong density gradients close to the boundary of the 
initial confinement. Since the latter ceases to act at t — 0, this leaves at t = + 
instantaneously a region R\ < r < R2 that is devoid of particles but in which the new 
external potential is essentially zero. This leads to a collective diffusion of the particles 
towards the boundaries set by the new potential. Correspondingly, the high density 
peaks decrease rapidly and leak outward. In the inner region, r « of the profile, the 
dynamics is much slower and the relaxation to the final plateau there takes place at the 
end of the process, causing thereby the final development of the new, weaker correlation 
peaks close to the location of the boundary, r < R 2 . 

The compression process, depicted in Fig. Efb) runs very differently. There, the 
initial 'closing' of the potential from Ri to R2 < R± leaves at t = + a region of 
high density at R\ < r < R2 that now finds itself within a strongly repulsive external 
field. There is an extremely rapid shrinking there, accompanied by the development of 
very high correlation peaks that actually 'overshoot' in height with respect to the final 
equilibrium profile. Initially, the region in the centre of the sphere remains unaffected 
and only as the high peaks start diffusing does material flow toward the centre and at 
the latest stage of the dynamics the profile at r = reaches its new equilibrium value. 

In order to quantify better this asymmetry and also extract characteristic time 
scales for the two dynamical processes, we consider the second moment of the density, 
m2(t), defined through 



The quantity m 2 (t) is a quantitative measure of the spread of p(r,t) around the centre 
of the external field and its time evolution is shown in Fig. El Let the superscript 




(12) 
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'+' denote the expansion and the superscript '— ' the compression process. Obviously, 
it holds m 2 (0) = ra^oo). We notice that for both processes m 2 (t) is a monotonic 
function of t but some interesting differences arise when one fits the two curves by 
analytic functions, shown as lines in Fig. El The expansion can be very accurately 
described by a single exponential: 

m+(t) = m+(0) + [m+(oo) -m+(0)][l - exp(-t/r + )], (13) 

with the characteristic time scale r + = 0.287 Tb- However, a double-exponential fit is 
necessary to parameterise the compression process, namely 

mj(i) = mj (oo) + A" exp(— i/rf ) + [m^"(0) — mj(oo) — A - ] exp(— t/r-f), (14) 

with the fit parameter A~ = 0.240 <r~ 2 and the two characteristic time scales r-f = 
0.036 tb and r 2 ~ = 0.189 tb- Since rf 2 < r+ ) it follows that the compression process is 
at any rate faster than the expansion one. The occurrence of the two distinct time scales 
r-f <C r 2 ~ in the compression requires some explanation. The fast process that takes 
place at times t ~ corresponds to the abrupt shrinking of the profile at the wings 
of the distribution and is caused exclusively by the 'closing' of the external field. This 
is the same mechanism that brings about the overshooting of the density peaks. Once 
this is over, diffusion within the now already confined system takes place and the second 
characteristic time scale, r 2 ~, is solely determined by the interaction potential V(r) and 
the average particle density. In the expansion process, the first mechanism is absent thus 
a single time scale, r + , shows up, which is of intrinsic origin exclusively. Since stronger 
density gradients occur during the compression than during the expansion process, even 
the larger of the two time scales of the compression, r 2 ~, is smaller than r + . The denser 
the system is, the faster the collective diffusion towards equilibrium. 

Finally, we turn our attention to the slab confinement. The results from theory 
and simulation are shown in Fig. 0] Once more it can be seen that the DDFT offers 
an excellent description of the dynamics of the system. The same asymmetry between 
expansion and compression that was seen in the spherical confinement also shows up for 
the case of the slab, including the overshooting of the peaks during the compression 
process. In addition, the density profiles develop during their evolution secondary 
oscillations that are also very well reproduced by the theory. 

In summarising, we have demonstrated that the dynamical density functional theory 
of Marini Bettolo Marconi and Tarazona [H |H], when supplemented by an accurate 
equilibrium density functional, can provide an excellent description of out-of-equilibrium 
dynamics of colloidal systems at the Brownian time scale. The accuracy of the DDFT 
formalism has already been successfully tested for the system of one-dimensional hard 
rods [4 , for which the exact density functional F[p] is known. To the best of our 
knowledge, this is the first study of the validity of DDFT in three dimensions. As 
the phenomenology in 3d is much richer than in Id, including the possibility of phase 
transitions, many intersecting ways for future applications open up. 
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Figure 4. DDFT (solid lines) and BD (noisy lines) results for the time development of 
the linear density profile p(z) in a slab confining potential V^(z, t) with (a) Z\ = 4.0 
and Zi = 5.0 and (b) Z\ — 5.0 and Zi = 4.0. The shown profiles are for the times 
t = Q,tx = 0.06, t 2 = 0.24 and t 3 = 2.0 (in units of r B ). 



Acknowledgments 

Discussions with Andrew Archer, Bob Evans, Wolfgang Dieterich and Pedro Tarazona 
are gratefully acknowledged. This work has been supported by the Deutsche 
Forschungsgemeinschaft through the SFB TR6. 



References 



[1] Evans R 1979 Adv. Phys. 28 143 

[2] Likos C N 2001 Phys. Rep. 348 267 

[3] Dhont J K G 1996 An Introduction to Dynamics of Colloids (Amsterdam:Elsevier) 

[4] Marini Bettolo Marconi U and Tarazona P 1999 J. Chem. Phys. 110 8032 

[5] Marini Bettolo Marconi U and Tarazona P 2000 J. Phys.: Condens. Matter 12 A413 

[6] Dieterich W, Frisch H L and Majhofer A 1990 Z. Phys. B 78 317 

[7] Lang A, Likos C N, Watzlawek M and Lowen H 2000 J. Phys.: Condens. Matter 12 5087 



Mean-field dynamical density functional theory 



9 



[8] Likos C N, Lang A, Watzlawek M and Lowen H 2001 Phys. Rev. E 63 031206 

[9] Louis A A, Bolhius P G and Hansen J P 2000 Phys. Rev. E 62 7961 

[10] Archer A J and Evans R 2001 Phys. Rev. E 64 041501 

[11] Archer A J and Evans R 2002 J. Phys.: Condens. Matter 14 1131 

[12] Archer A J, Likos C N and Evans R 2002 J. Phys.: Condens. Matter 14 12031 

[13] Archer A J, Evans R and Roth R 2002 Europhys. Lett. 59 526 

[14] Louis A A, Bolhuis P G, Hansen J P and Meijer E J 2000 Phys. Rev. Lett. 85 2522 

[15] Allen M P and Tildesley T J 1989 Computer Simulation of Liquids (Clarendon Press: Oxford) 

[16] Ermak D L 1975 J. Chem. Phys. 62 4189 



